{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "from osgeo import ogr, osr\n",
    "\n",
    "from MalardClient.MalardClient import MalardClient\n",
    "from MalardClient.DataSet import DataSet\n",
    "\n",
    "client = MalardClient(\"MALARD-PROD\")\n",
    "\n",
    "\n",
    "proj4 = client.getProjection(DataSet( \"test\", \"swath_c_GrIS_OIB\", \"greenland\"  )).proj4\n",
    "print(proj4)\n",
    "\n",
    "\n",
    "daShapefile = \"/data/eagle/project-cryo-tempo/data/masks/greenland/NoPeripheryMerged/GRE_Basins_IMBIE2_v1.3_NoPeriph.shp\"\n",
    "#daShapefile = \"/data/puma1/scratch/mtngla/chucach/Chucach-polygon.shp\"\n",
    "\n",
    "#kmz file\n",
    "#daShapeFile = \"/data/mouse1/team/jon/Anarctica/Coastline_Antarctica_v02.shp\"\n",
    "\n",
    "\n",
    "# Save extent to a new Shapefile\n",
    "outShapefile = \"icesheet_noperiph.shp\"\n",
    "\n",
    "driver = ogr.GetDriverByName('ESRI Shapefile')\n",
    "\n",
    "dataSource = driver.Open(daShapefile) # 0 means read-only. 1 means writeable.\n",
    "feature = None\n",
    "\n",
    "layer = dataSource.GetLayer(0)\n",
    "\n",
    "for i in range(0, layer.GetFeatureCount()):\n",
    "    f = layer.GetFeature(i)\n",
    "    \n",
    "    print(f.GetField(0))\n",
    "#    if str(f.GetField(1)).startswith(\"LRM\"):\n",
    "#        print( f.GetField(1))\n",
    "    \n",
    "#    if f.GetField(1) == \"SARIN over Greenland Coast\":    #SARIN over Greenland Coast, LRM over Antarctica ice sheet, LRM over Greenland ice sheet,  \n",
    "#        feature = f\n",
    "#        print(f.GetField(1))\n",
    "#    if f.GetField(1) == \"SARIN over Greenland Coast\":\n",
    "#        print(f.GetGeometryRef().GetEnvelope())\n",
    "\n",
    "\n",
    "#greenlandExtent = (-81.265, -6.8349, 59.1676, 84.2163)\n",
    "\n",
    "        \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "sourceprj = layer.GetSpatialRef()\n",
    "\n",
    "srs = osr.SpatialReference()\n",
    "srs.ImportFromProj4(proj4)\n",
    "\n",
    "targetprj = srs\n",
    "transform = osr.CoordinateTransformation(sourceprj, targetprj)\n",
    "        \n",
    "        \n",
    "outDriver = ogr.GetDriverByName(\"ESRI Shapefile\")\n",
    "\n",
    "# Remove output shapefile if it already exists\n",
    "if os.path.exists(outShapefile):\n",
    "    outDriver.DeleteDataSource(outShapefile)\n",
    "\n",
    "# Create the output shapefile\n",
    "outDataSource = outDriver.CreateDataSource(outShapefile)\n",
    "outlayer = outDataSource.CreateLayer(\"icesheet_noperiph\", targetprj,geom_type=ogr.wkbPolygon)\n",
    "\n",
    "# Add an ID field\n",
    "idField = ogr.FieldDefn(\"id\", ogr.OFTInteger)\n",
    "outlayer.CreateField(idField)\n",
    "\n",
    "# Create the feature and set values\n",
    "for i in range(0, layer.GetFeatureCount()):\n",
    "    feature = layer.GetFeature(i)\n",
    "\n",
    "    transformed = feature.GetGeometryRef()\n",
    "    transformed.Transform(transform)\n",
    "    \n",
    "    def replacePoint(wkt, x, y ):\n",
    "        strToReplace = \"{} {}\".format(x, y)\n",
    "        strReplacement = \"{} {}\".format(x, y + 100)\n",
    "        \n",
    "        return wkt.replace(strToReplace, strReplacement)\n",
    "       \n",
    "    gb = transformed.Buffer(0)\n",
    "    \n",
    "    geom = ogr.CreateGeometryFromWkb(gb.ExportToWkb()) \n",
    "    \n",
    "    print(geom.IsValid())\n",
    "    \n",
    "    defn = outlayer.GetLayerDefn()\n",
    "    feat = ogr.Feature(defn)\n",
    "    feat.SetGeometry(geom)\n",
    "    outlayer.CreateFeature(feat)\n",
    "    \n",
    "\n",
    "# Save and close DataSource\n",
    "inDataSource = None\n",
    "outDataSource = None\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "#Warning 1: Self-intersection at or near point 182103.98757600426 -2898036.0123973158\n",
    "#Warning 1: Self-intersection at or near point 540934.02211438434 -1894356.1721276832\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
